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Abstract 

We solve the 2D Dirac equation describing graphene in the presence of a hnear vector potential. 
The discretization of the transverse momentum due to the infinite mass boundary condition reduced 
our 2D Dirac equation to an effective massive ID Dirac equation with an effective mass equal to the 
quantized transverse momentum. We use both a numerical Poincare Map approach, based on space 
discretization of the original Dirac equation, and direct analytical method. These two approaches 
have been used to study tunneling phenomena through a biased graphene strip. The numerical 
results generated by the Poincare Map are in complete agreement with the analytical results. 
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1 Introduction 



Graphene, a single layer of carbon atoms laid out in a honeycomb lattice, is one of the most interesting 
electronic systems discovered in recent years [H |2] . It differs from conventional two dimensional elec- 
tron gas (2DEG) systems in that the low energy physics is governed by a massless Dirac Hamiltonian 
rather than the more common form used for semiconductors, characterized by an effective mass and 
a band gap. 

The tunneling of Dirac fermions in graphene has already been verified experimentally [3] , which in 
turn has spurred an extraordinary amount of interest in the investigation of the electronic transport 
properties in graphene based quantum wells, barriers, p-n junctions, transistors, quantum dots, su- 
perlattices, • • • etc. The electrostatic barriers in graphene can be generated in various ways [H E], by 
applying a gate voltage, cutting it into finite width nanoribbons and using doping or otherwise. On 
the other hand, magnetic barrier could in principle [H [6l [7] be realized with the creation of magnetic 
dots. In the case of graphene, results of the transmission coefficient and the tunneling conductance 
were already reported for the electrostatic barriers, [H O [U [71 [8] , magnetic barriers, [H [9l [10] and 
potential barrier [TTl [T2l [T3] . 

The fact that in an ideal graphene sheet the carriers are massless gives rise to Klein paradox, which 
allows particles to tunnel through any electrostatic potential barriers, that is the wavefunction has an 
oscillatory tail outside the electrostatic barrier region. Hence this property excludes the possibility to 
confine electrons using electrostatic gates, as in usual semiconductors. Thus to enable the fabrication 
of confined structures, such as quantum dots, we need to use other type of potential coupling such as 
the vector potential coupling considered in the present work. More precisely, we will use two methods, 
one numerical (Poincare map) [13] and the other analytical, in the prospect of studying the tunneling 
of Dirac fermions through a biased graphene strip. 

The paper is organized as follows. In section 2, we describe our theoretical model Hamiltonian 
and apply the Poincare map approach, based on the space discretization of the effective ID Dirac 
equation. In section 3, we expose the direct analytical approach to solve the same problem. In section 
4, we proceed to discuss the numerical implementation of our approaches to a specific model potential, 
the linear potential which generates a static electric field, and make a comparative study between the 
two approaches. 

2 Poincare map 

Before we embark on the two approaches mentioned above, we would like to describe mathematically 
our system of massless Dirac fermions within a strip of graphene characterized by a very large length 
scale, and a width W in the presence of the applied linear potential V{x) between x = and x = L. 
So our system is composed of three major regions: the extremes (I) and (III) contain intrinsic graphene 
free of any external potentials and an intermediate region (II) subject to the applied linear potential 
V{x). Graphene band structure has two Fermi points, each with a two-fold band degeneracy, and can 
be described by a tight binding Hamiltonian describing two interlacing honeycomb sublattices. At 
low energies this Hamiltonian can be can be described by a continuum approximation to the original 
tight binding model which reduces to the two dimensional Dirac equation with a four-component 
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envelope wavefunction whose components are labeled by a Fermi-point pseudospin = ±1. Specifically, 
the Hamiltonian for one-pseudospin component for the present system can be written as 

H = VF^ -p + Vix) (1) 

where vp ~ 9.84 x lO^m/s is the Fermi velocity and B = {(Tx,(Ty) are the Pauli matrices. Hereafter we 
set our units such that vp = h = 1. The linear potential V{x) has the following form 

-Fx + Vo, 0<x<L 
0, otherwise 

where F = ^ is the strength of the static electric field. This potential configuration is shown in 
Figure 1 below. 





Figure 1: Discretization of the linear potential V{x). 

Our system is supposed to have a finite width W with infinite mass boundary conditions for the 
wavefunction at the boundaries y = and y = W along the y-direction |1H I15j . This boundary 
condition results in a quantization of the transverse momentum along the y-direction, which gives 



k - — 



/ = 0,1,2 • 



(3) 



One can therefore assume a spinor solution of the following form 'tp^[x,y) = ( (/>2(x) I e'^^y^ 



t 



where the superscript j = 1,11,111, indicates the space region while the subscripts indicate the two 
spinor components. Thus our problem reduces to an effective ID problem whose Dirac equation can 
be written as 



V{x) - e 
-1 + ky V{x 



( 4>{{x) \ 



0. 



(4) 



" ^ \ -#2(^) / 

Due to the space dependence of the potential V{x) we make the following transformation on our spinor 



components to enable us to obtain Schrodinger like equations for each component, Xi 



2 



and xi ~ ^ (jPi ~ "^Cj ' which obey the coupled stationary equations. These are 



dx 



± i {V{x) - e) x{o{^) =F ikyxiiix) = 0. 



(5) 



Each spinor component xi 2 ^e shown to satisfy the following uncoupled second order differential 
equation 



d 



^^2X1,2 (^) + (^) + (^) - - ) xi,2 (^) = 0. 



(6) 



In this section we will apply the Poincare map approach to solve the above effective ID Dirac 
equation. In this approach we start by subdividing the potential interval L into + 1 regions (Figure 
1). In every n-th region we approximate the linear potential by a constant value Vn = V{xn) where 
Xn = nh and h = j^rj- Hence, the Dirac equation in each region (n), defined by h{n — 1) < a; < hn, 
can be easily solved for the piece-wise constant potential. For simplicity, we chose the incident wave 
propagating from right to left and apply the continuity of the spinor wavefunctions at the boundary 
separating adjacent regions. The general solutions of equation ([6|) in the n-th region where V{x) = Vn 
are given by 

' 1 \ ^ „ / 1 ' 



(7) 



with kn = — Vn)'^ — k1, the complex number Zn is defined by = ^ = sgn (e — Vn) rr^ — -7. 
In order to obtain the relationship between ipn+i and ipn we apply continuity of at the boundary 
X = Xn (Figure 2). This leads to 



Mn{Xn) y ^ j =Mn+l{Xn) 

Also Mn+i{xn+i) and Mn+i{xn) are related by 

Mn+l{Xn+l) = Mn+l{Xn)Sn+l-, Sn+1 = 



Bn+1 



(8) 



-ihk„+i 





(9) 
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Figure 2: Solutions of the ID Dirac equation in two consecutive regions, continuity of spinors is applied 
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Using the above results we can write the desired Poincare map as 

1pn+l{Xn+l) = Tn-fpniXn) 



(10) 



where we have defined a simpHfied notation by ipn = '4'n{xn) and r„ = M„_|_i(x„+i)5„+iMjj_,^^(a;„+i), 
or more exphcitly 



^(i+„) + ^(i+n) V -e-^'=(i+-) + e-^'=(i+-) 4+„)e-^^(i+") + Z(i+„)e-*'^(i+" 



(11) 



To make use of the above Poincare map in solving our scattering problem we need to define our 
incident, reflected and transmitted waves. For 2; < where V = Q (region I), we can use for our 
transmitted spinor evaluated at n = 0, the suitably normalized form 

*° - ( ) ■ (-) 

This is juste the value of the transmitted wave at the zeroth site, x = (n = 0), which is given by 

i^L = Ao(^ (13) 

On the other side, for x > h{N + 1) where V = Q (region III), we have both incident and reflected 
spinor waves. Just outside the potential region on the right hand side in the {N + 2)-th region the 
spinor wave can be written as 

i^R = ^ e-*'^"^ + i?jv+2 ^ (14) 

Hence to evaluate the transmission amplitude all we need is to find using the above recursive 

scheme. Our strategy now is to express Apfj^2 in terms of V'Af+i and '0Af+2) the two end point spinors. 
This can be easily done using our previous relationships and leads to 

^^+2 = ^(^_^2.Mo) (1 ) (^^+2 - e^"'°^I^N+l) . (15) 

From the above notation we can easily define the transmission amplitude as follows 

An+2 ^ ^ 

Summing up, our numerical procedure requires first that we iterate the Poincare map pO|) to obtain 
the end point spinors, ^tv+i and in terms of the normalized transmitted spinor. These spinors 

will then be injected in (jlSp and ()16p to determine the transmission amplitude. The transmission 
coefficient is given by T = |i| . The numerical implementation of this scheme in the case of linear 
vector potential will be done in section 4. 

3 Analytical method 

Let us now solve analytically the effective ID Dirac equation or equivalently equation ([6]) in the 
presence of a quantum barrier (region II). Our objective is to find the transmission coefficient for a 
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Dirac fermion scattered by a linear potential and then compare our results with those found in previous 
section using the Poincare map method. The solution of equation ([6]) in region I and 1 1 1 are given by 

where r and t are the reflection and transmission amplitudes, respectively. The wave vector kx = 



£^ — and the complex number z is defined by z = sgn(e)(A:2: + iky)/yjk'^ + ^y- region II the 
general solution can be expressed in terms of the parabolic cylinder function [161 [17] as 

x'lix) = aD,., (^^e'-/\Fx + E)j + f3D_, (^-^e-'^^/^Fx + E)j (18) 
where u = E = e — Vq, a and /5 are constants. Substituting (fT8]l in ([5]) gives the other component 



X^^{x) = ^2{E + Fx)D^^ (^-^|.e-'^/^(Fx + E)j + ^/2Fe*^/4L»_^+i [- ^ ^e-^^'^{Fx + E) 

-^V2Fe~^^/^Dy (\^e^^/^{Fx + E)) (19) 



The components of the spinor solution of the Dirac equation ([T]) in region II can be obtained from (|18|) 
and (|19p where 4>^i{x) = y^l + iy^2 ^'^id '/'2(^) ~ ~ ^X2- This results in 

where the function a^(x) and 6^(x) are given by 

a±(x) = Z),_i ( y|e-/4(Fx + E)\ T ^e^/^Z?, (Me'^/\Fx + i?) ] 



6±(x) 



±^{-2iE ±ky- 2iFx)D_^ l-^e-'^/^{Fx + EU . (21) 



The coefficients r, a, (3 and t are determined from the continuity of the spinor wavefunctions at the 
boundaries x = 0,L, that is ^/^'(x = 0) = t(j^^{x = 0) and Tp^\x = L) = -(/^'"(x = L). The transmission 
coefficient through the linear potential is obtained from T = \t\ where the corresponding amplitude 
t is obtained from the aforementioned boundary conditions. It is given by 

e-ik.L ^ [b+{L)a~{L) - b'{L)a+{L)\ 



[6+(0) + z6-(0)] [a-{L) - za+{L)] - [a+(0) + za-(O)] [b-{L) - z6+(L)] ' 

4 Results and discussion 



(22) 



In this section we implement our previous Poincare map and analytical approaches to a nanoribbon 
system subject to an electric potential of strength Vq = 10, 20 and a field region of length L = 3, 10 so 
that the resulting static electric field strength is given hy F = Vq/L = 10/3, 2, respectively. In Figure 
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3 we show the transmission as a function of energy for a transverse momentum ky = 1. The sohd 
hnes corresponds to the exact transmission derived in section 3 and given by equation (j22p while the 
dashed hnes are generated by our Poincare map for N = 200 iterations, the agreement is just perfect. 




10 15 20 
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Figure 3:Transmission coefficients T versus energy efor L = 3, 10, Vq = 10,20 and ky = 1. 

Figure 3 shows the concordance between the results generated by the analytical and Poincare map 
method we adapted. We note that below a certain critical energy e = ky the transmission is almost 
zero, then it starts oscillations whose frequency increases with L, the size of the region subject to the 
electric field. The transmission increases with L and reaches unity for energies above Vq + 2ky. 
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Figure 4:Transmission coefficients T versus Vq for L = 3, 10, e = 10, 20 and ky = 1 

Figure 4 shows the transmission as a function of the strength of the applied voltage, total transmission 
is observed for small values of Vq less than the energy of the incident fermion. It then decreases sharply 
for Vq > e — 2ky until it reaches a relative minimum and then begins to increase in an oscillatory 
manner. We notice in both Figures 3 and 4 that the amplitude of oscillations and period increase as 
we decrease the size of the electric field region, L. 
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Figure 5:Transmission coefficients T versus energy e for L = 3, Vq = 20 and different values of ky. 

Figure 5 shows that the effect of the transverse momentum ky on transmission, is antagonistic to that 
of length L. But it should be pointed out that the number of oscillations increases as ky decreases 
and the curves for different values of ky do not intersect. 

Now we would point out that our effective ID massless Dirac equation is equivalent to a massive 
one with an effective mass equal to the transverse quantized wave vector ky. For this purpose we would 
like to consider a unitary transformation, which enable us to map the effective ID (equation (jH)) into 
a ID massive Dirac equation. Such a unitary transformation does not affect the energy spectrum or 
the physics of the problem. We choose a rotation by 7r/4 about the y-axis, U = e ^'^'^ . Thus, the 
transformed Hamiltonian and wavefunction read 

= 0, 4,,(x) = i7^,-,_,(x). (23) 

V i'iix) j 

which is identical to a ID massive Dirac equation with an effective mass m = ky. To check the validity 
of this assertion numerically we show in Figure 6 the transmission as a function of energy as generated 
by the exact analytical result (j22p . the Poincare map (jl6p and the ID massive Dirac equation with 
an effective mass m* = ky in ()23p . We see from this figure that the three curves coincide to the point 
that we cannot even distinguish between them. This lead us to include an inset in Figure 6 showing 
each figure translated for ease of comparison purposes. 

1.0 
0.8 
0.6 
0.4 
0.2 
0.0 

10 20 30 40 50 60 

Figure 6 :Transmission coefficients T{e) for L = 1, Vq = 40 and ky = m* = 

This last figure confirms, numerically, the equivalence between a one-dimensional system of Dirac 
fermions with mass and a two-dimensional system of massless Dirac fermions constrained along the 
y-direction by an infinite mass boundary condition, hence forming a graphene nanoribbon. The 
transverse component of the wave vector, ky, played the role of an effective mass [T5] in the resulting 
effective ID Dirac equation. 
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